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Abstract 

We derive a model-independent integral formula for chiral susceptibility and attempt to present 
a continuum model study of it within the framework of Dyson-Schwinger Equations. An appro- 
priate regularization is implemented to remove the temperature-independent quadratic divergence 
inherent in this quantity While it demonstrates a second-order phase transition characteristic in 
the chiral limit, the result obtained supports a crossover at physical current quark masses, which 
is in good agreement with recent lattice studies. 
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I. INTRODUCTION 



Quantum Chromodynamics (QCD) underlies our current understanding of strong interac- 
tion and explains why and how the fundamental constituent, namely quarks and gluons are 
bound into nucleons, pions, etc. The QCD vacuum is rather complex and intriguing: while 
it features the overwhelming spontaneous chiral symmetry breaking and color confinement, 
it is believed to be heated and undergo a transition at high enough temperature 
As a result, quarks and gluons are no longer confined inside hadrons and a new state of 
strongly interacting matter, the so-called quark-gluon plasma (QGP) is expected to form. 
Meanwhile, the chiral symmetry gets restored. It is a goal to detect such a new state of 
matter for the on-going heavy-ion collision experiments at the BNL Relativistic Heavy-Ion 
Collider (RHIC) and the near future Large Hadron Collider (LHC). It is widely speculated 
that this QCD-related transition may also be of direct relevance for the evolution of the 
early Universe l|, |5| 

Although lattice QCD has established the existence of such a thermal transition, the 
nature of this transition (its order, or, whether is is a real phase transition) has been a long 
time debate [l]. It is not until recently that Y. Aoki and his collaborators, by exploiting 
quite demanding and powerful lattice simulations, assert unambiguously that QCD evolves 
smoothly with temperature; that is, there is no thermodynamic phase transition for three 
flavors (u, d, s-quark) with physical current masses, but instead a smooth crossover 
Such a conclusion rules out some scenarios and speculations over the evolution of the early 
Universe jl|. In fact, the nature of the QCD thermal transition would lay an imprint on 
the Universe's evolution and thus affect our understanding of the latter. For example, the 
previously prevailing strong first-order transition scenario predicts inhomogeneous nucle- 
osynthesis and later formation of cold dark matter clumps Nowadays, people tend 
to believe that the nature of QCD thermal transition depends on the quark content and 
current quark masses . In the chiral limit, the chiral phase transition is likely to be 
of second order for two flavors and first order for three flavors. For infinite or very large 
current quark masses, the transition turns out to be a first-order deconfinement transition. 
For intermediate current masses, one may end up with a crossover, instead of a real phase 
transition, meaning that thermodynamic quantities would show a rapid change, as opposed 
to a jump, in a narrow temperature range. 



2 



As is well known, first principle lattice simulations provide the prior choice to address 



its intrinsic nonperturbative nature. In 
6] achieved a large scale computing and 



the thermal QCD transition problem, because o: 
this regard, recently Y. Aoki and collaborators 
identified the nature of three-flavor QCD finite temperature transition as a smooth crossover 
at physical quark masses. On the other hand, other than lattice simulations, continuum 
nonperturbative model study can be a complementary approach to thermal QCD transition. 



Following Y. Aoki et al. in Refs. j^, [6|, we present a continuum model study of the 
interesting chiral susceptibility in the present work, from whose behavior as temperature 
varies we identify the nature of two-flavor QCD thermal transition in the case of physical 
current quark masses as well as in the chiral limit. 

We first derive a model-independent analytical formula, which expresses chiral suscepti- 
bility as a integral in terms of dressed propagators and dressed scalar vertex. The quadratic 
divergence inherent in this susceptibility is then demonstrated and an appropriate renormal- 
ization procedure is specified. In the section that follows, this result is then illustrated by 
a model calculation based on an effective interaction in the framework of rainbow Dyson- 
Schwinger equation (DSE) and ladder Bethe-Salpeter equation (BSE) approach. The no- 
torious quadratic ultraviolet divergence is appropriately subtracted by a delicate numerical 
procedure. The results obtained are carefully analyzed and a crossover and a second-order 
phase transition are respectively identified for the cases of physical and zero current quark 
masses. The final section is devoted to the conclusion and discussion. Throughout this 
paper, we work with Euclidean space metric: {7^,7^} = 25^ u . 



II. DERIVATION AND RENORMALIZATION 



A. Derivation of an integral formula for chiral susceptibility 



Chiral susceptibi 



ity has interested quite a few 



much investigation 0, y, y, y, q, q, q, q. 
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m 



attice QCD practitioners and attracted 
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13, 



201 ] . It is defined as 



T d 2 d - 

x = — nz? lo g z = (-)q— < H' >> 



V dm 2 



dm 



(1) 
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where Z denotes the QCD partition function at finite temperature and vanishing chemical 
potential. 

Z = J V^VAexp{- J dr J d 3 s#/(7 • d + m f + ig^ ■ A a )^ f + \f^F^\}. (2) 

Here, vp, ip and are subject to anti-periodic and periodic conditions, respectively; and 
we leave the gauge fixing term and the ghost field term to be understood. Note that chiral 
susceptibility is defined for two light flavors and hence m in Eq. (1) is supposed to be the 
degenerate current mass of u, d quarks. The second line of that equation also indicates that 
the chiral susceptibility measures the response of chiral condensate (the order parameter) to 
a small perturbation of the parameter responsible for explicit breaking of chiral symmetry 
(the current quark mass). 

Formally, using quark propagator we can always write the quark condensate as 

^- 4 tr y G(p,m), (3) 

where tr 1 denotes trace over Dirac indices of the quark propagator, N c = 3 is the color 
factor and Nf — 2 denotes two degenerate light flavors. Substituting Eq. (3) into Eq. (1) 
and adopting the identity 

^) = - G (p, m) ^^) G (p, m ), (4) 

we have 

X = (~)N c N f J ^tr^G(p,m) dG ~^ m) G{p,m)]. (5) 

We consider the mass m as a constant background scalar field coupled to the quark fields by 
the term mipip. Then G(p, m) is the dressed quark propagator in the presence of such a back- 
ground field and the derivative of its inverse with respect to the current quark mass yields 
the so-called one-particle- irreducible (1PI) dressed scalar vertex (a three-point correlation 
function) 

dG-\p,m) 

— = r(p,0;m), (6) 

where p is the relative momentum. This relation may be called "scalar Ward identity" . Note 
that the total momentum of the dressed scalar vertex vanishes because the background scalar 
field m is a constant that takes no momentum. Substituting Eq. (6) into Eq. (5) gives 

—Ltr^Gfa m)T(p, 0; m)G(p, m)\. (7) 



This is the integral formula for chiral susceptibility at zero temperature. Diagrammatically, 
the right- hand- side of Eq. (7) expresses the Feynman diagram for the color-singlet scalar 
vacuum polarization at zero total momentum. This clearly shows that chiral susceptibil- 
ity at zero temperature equals the color-singlet scalar vacuum polarization at zero total 
momentum. 

According to finite-temperature field theory (see, e.g., 181]). the corresponding finite 
temperature version is obtained by replacing the integration over the fourth component of 
momentum with summation over fermion Matsubara frequencies: 

+00 



X(T) = (- 



)N c N f T j^tr 1 [G( Pn ,m)T(p n ,0;m)G( Pn ,m)l (8) 



where p n = (p,u> n ) with Matsubara frequencies uj n = (2n + 1)ttT. Here it should be noted 
that in the chiral limit the above expression for the chiral susceptibility reduces to the 
corresponding expression (Eq. (7)) in Ref. [jij. 

So we have derived a model-independent integral formula which expresses the chiral sus- 
ceptibility in terms of the dressed quark propagator and the dressed scalar vertex, both of 
the latter being basic quantities in quantum field theory (QFT). The DSE-BSE approach 
provides a desirable framework to calculate these quantities non-perturbatively and hence 
the chiral susceptibility. 



B. Quadratic divergence and renormalization 

There resides an additive quadratic ultraviolet divergence in the chiral susceptibility, 
which can manifest itself even in the chiral limit 0,0. To see this unambiguously, we 
calculate the chiral susceptibility in the case of free quark gas, because all quantities tend to 
the free case counterparts in the ultraviolet limit (asymptotic freedom) and all divergences 
originate from the ultraviolet limit. In this case, the dressed scalar vertex reduces to the 
bare one, i.e. 

r(p n ,o)-i, (9) 

and the free quark propagator in the chiral limit reads 

G free ( Pn ) = — !— . (10) 

11 -Pn 
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Substituting Eqs. (9) and (10) into Eq. (8), one gets 

d 3 p 



X 



/-(r) = 4 W E/(^^ (ID 

n=— oo v ' 

With the aid of the identity 



+oo 1 
X 1 



CO 2 + X 2 2 
n=— oo " 



where np{x) = eX /r +l is the Fermi-Dirac statistics function, we have 

X*~(T) = 2N c N f J - 2n F (\p\)). (13) 

The second part of the integrand can be integrated out and gives a purely temperature- 
dependent finite result: ^ T 2 . On the contrary, the first part yields an additive quadratic 

divergence. This is in striking contrast to the free quark-number susceptibility, which is a 



definitely finite quantity proportional to T 2 and therefore vanishes at zero temperature 19]. 

Actually, there always contains an inherent additive quadratic divergence in chiral suscep- 
tibility. One dose not really have to address the issue of the regularization of this quadratic 
divergence, when one works with a hard cut-off in the numerical integration (e.g., lattice 
simulations), and is concerned only with the temperature-dependent behavior of the suscep- 



tibility, not so much on its absolute value 2jJ . However, a correct renormalization procedure 



is important and indispensable if one wants to get correct physics in the continuum limit 
fl. In Ref. ljj, we identified the disconnected part of the chiral susceptibility and argued 
that it is of major relevance and free from quadratic ultraviolet divergence. However, we 
were unable to derive the disconnected chiral susceptibility directly from a first-principle 
approach such as the functional integral method. Looked at from this point, the identifi- 
cation of the disconnected chiral susceptibility is not yet absolutely "waterproof" and solid 
enough. So in this paper we follow the renormalization procedure in Ref. which is to 
subtract the zero temperature susceptibility and study the difference between T ^ and 
T = 0. That is to say, we define the renormalized chiral susceptibility by 

Xr{T)=x(T)-X (14) 

where x(T) and x are provided by Eq. (8) and Eq. (7), respectively. The justification for this 
renormalization is two-fold: firstly, the quadratic divergence is additive and temperature- 
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independent and hence can be really removed by subtracting the zero-temperature suscep- 
tibility; secondly, this subtraction will never affect the genuine temperature effects of chiral 
susceptibility that we are interested in. 

Having specified the renormalization procedure, in the next section, we will focus on 
the model numerical calculation of the renormalized chiral susceptibility Xr{T) within the 
framework of DSE-BSE. 



III. MODEL NUMERICAL CALCULATION 

Dyson- Schwinger equations provide a non-perturbative, continuum approach for the ex- 
ploration of strong interaction physics [2lJ. Derived from QCD's Euclidean space generating 
functional, they constitute a tower of enumerable coupled integral equations whose solutions 
are the n-point Green's functions. A consistent truncation scheme that makes the DSE cal- 
culations tractable is the so-called rainbow-ladder approximation. Such a rainbow-DSE and 
ladder-BSE approach has been used extensively in the investigation of spontaneously chiral 



symmetry breaking and confinement 22j and also found successful applications in calcu- 



lating light pseudo-scalar and vector meson observable s [231 . [24| as well as in describing 
strong interaction at finite temperature and/or density [25j. The subsequent model calcu- 
lation of chiral susceptibility will also be conducted systematically and consistently in the 
rainbow-DSE and ladder-BSE framework with an effective gluon propagator. 



A. The zero-temperature chiral susceptibility 

We first calculate the zero-temperature chiral susceptibility with Eq. (7). This involves 
the calculation of dressed quark propagator and dressed scalar vertex. The former is solved 
from the rainbow gap equation 

4 f d 4 q 

G-\p,m) =i~?-p + m + - J j^y9 2 D e JJ {p - gh M G(g, (15) 

and the latter satisfies the ladder inhomogeneous BSE (the total momentum of the dressed 
scalar vertex has been set equal to zero) 

r(p, 0; m) = 1 - ~ J AJLfD'J/ip - q)^G(q, m)T(q, 0; m)G(q, m) 7 „, (16) 
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where D e JJ {p—q) is the effective gluon propagator, which is usually a phenomenological input 
in practice (it should also be noted that the analytical structure of the gluon propagator has 
been explored from numerical solutions of coupled DSEs of quarks, gluons and ghosts and 
from fits to lattice data in recent literatures, see Ret. Q and references therein). 

At the moment, we pause for a while and give a brief proof that the rainbow-ladder 
truncation scheme of DSE-BSE respects the scalar Ward identity (Eq. (6)). Taking the 
derivative with respect to m on both sides of Eq. (15) and using again the identity Eq. (4), 
one arrives at 

= 1 - 3 J W jD " (f - MiM— ^Gh,m)7.. I") 

Comparing Eq. (16) and Eq. (17), one sees that r(p, 0;m) and dG~ 1 (p,m)/dm satisfy 
the same integral equation and therefore they must be equal, i.e., the scalar Ward identity 
holds at the level of the rainbow-ladder truncation. That is to say, the rainbow-ladder DSE- 
BSE model calculation respects the scalar Ward identity, whatever form the effective model 
gluon propagator takes. In fact, the rainbow-ladder truncation scheme of DSE-BSE itself is 
consistent in a broader sense. In addition to the scalar Ward identity, this truncation scheme 
preserves the Abelian vector Ward-Takahashi identity and the axial-vector Ward-Takahashi 
identity {23]. It is the preservation of the latter identity that guarantees an understanding 
of chiral symmetry and its dynamical breaking in the rainbow-ladder DSE-BSE framework 



without a fine-tuning of model-dependent parameters 



23] 



In this paper, we will employ the rank-2 confining separable model gluon propagator 

g 2 D*JJ(jp - q) - 5, u D(p 2 , q\ V - q) = 5^[D f (p 2 )f (q 2 ) + DJ^p ■ qftf)} (18) 

with fi(p 2 ) = exp(—p 2 /A 2 ),i = 1, 2. It is found to be very successful in describing light flavor 
(including u, d, s-quarks) pseudo-scalar and vector mesons with parameters A = 0.638 GeV, 



27 



Ai/A = 1.21, D A 2 = 260.0, D X K\ = 130.0 and m = 0.0053 GeV 

To proceed our model calculation, note that the quark propagator can be decomposed in 
terms of two independent Lorentz structures (no 75-related structure appears in pure strong 
interaction) 

G _1 (p, rn) = «7 • pA(p 2 , to) + B(p 2 , m) (19) 
Upon substituting Eqs. (18) and (19) into Eq. (15), the latter is converted into two coupled 
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integral equations (s = p 2 ) 

a(m) .AT ds s 2 hm + a(m)A(s)} 

1 ) 2Att 2 J s[l + a(m)/ 1 (s)] 2 + [m + 6(m)/o(s)] 2 1 J 

6fm) = — /°° rf s */o(*)[™ + ft(™)/o(*)] 

1 J Wo stl + aHA^p + fm + ftH/o^)] 2 1 J 

with A(p 2 ,m) = 1 + a(m)fi(p 2 ) and B(p 2 ,m) = m + b(m)fo(p 2 ). 

As for the dressed scalar vertex, its general form from Lorentz structure analysis reads 

r(p, 0; m) = F(p 2 , m) • 1 + 17 • pG(p 2 , m), (22) 



but only the first term provides the leading behavior [29j. We shall keep only this term 
for simplicity. Putting this ansatz V(p, 0;m) = F(p 2 ,m) ■ 1, the model gluon propagator 
Eq. (18) and the quark propagator Eq. (19) into Eq. (16), one has after completing some 
algebras 

r(p, 0; m) = F(p 2 , m) = 1 + a(m)/ (p 2 ), (23) 
where a(m) satisfies the following equation 

. . 16D f d A q r . 9 . , . \ » / ox-, q 2 A 2 (q 2 ,m) - B 2 (q 2 ,m) , . 

= — y <2^ /ofe f x [i + a(m)/ ° fe ,] x <24) 

We point out that the Gaussian factors fo{p 2 ) and fi(p 2 ) in the model gluon propagator 
provide sufficient ultraviolet suppression and justify an effective cut-off for momentum inte- 
gration in Eqs. (20,21,24). Therefore no multiplicative renormalization is needed. 

With the dressed quark propagator (A(p 2 ,m) and B(p 2 ,m)) and then the dressed scalar 
vertex (a(m)) solved, the zero-temperature chiral susceptibility (Eq. (7)) can now be written 
as 

x(m) . iNcN , J0. ll + aW/0(p2)] x w^J^ (25) 

and calculated out. Of course, x( m ) takes an additive quadratic divergence, which serves 
to cancel that one contained in the temperature-dependent chiral susceptibility. 

B. The finite-temperature chiral susceptibility 

In recent years, the rainbow-ladder truncation DSE models are extended to describe 



the finite temperature properties of QCD and its chiral phase transition 25j. As for the 



foregoing confining separable model, the extension to finite temperature is systematically 



accomplished by transcription of the Euclidean quark four-momentum via q — > q n = (q, oj t , 



where u n = (2n 
are introduced 
modeled as 
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- 1)ttT are the fermion Matsubara frequencies, and no new parameters 
. This means that the effective gluon propagator at finite temperature is 



9 2 D% f {p k - q n ) = S^DoMM) + DtMpDp, ■ Qnfi(ql)}. (26) 

with pi = ~p 2 + uj\. The subsequent calculations involving temperature-dependent dressed 
quark propagator and scalar vertex will be parallel to those of the zero-temperature case 
specified in the last subsection. 

Firstly, the finite temperature version of the rainbow-DSE for the dressed quark propa- 
gator Eq. (15) reads 

4 f d 3 p 

G 1 (p k , m) = «7 ■ p k + -T / —--g 2 D e JJ(p k - g n )7 M G(g n , m)^ v . (27) 

n=— oo ' 



Its solution now consists of three independent amplitudes due to the breaking of 0(4) sym- 
metry down to 0(3) symmetry |25j: 

G _1 (p„, m) = 27 • pA(p 2 n , m) + i>y 4 uj n C(p 2 n , m) + B(p 2 n , m). (28) 

Substituting Eq. (28) into Eq. (27), one easily obtains solutions of the form = 
1 + a(T,m)fi(p 2 n ), B(pl,m) = m + b(T,m)f (p 2 n ), and C(p 2 n ,m) = 1 + c(T,m)fi(pl,m), 
where a(T, m), b(T, m) and c(T, m) are functions of T satisfying the following three coupled 
nonlinear equations: 

a(T,m) = £ / 7^M)f[l + a^m)^)}^ 1 ^^) (29) 

n=— 00 ^ ' 

c(T,m) = ^T £ J ^J 1 (pl)u J 2 n [l + c(T,m)f 1 (p 2 n )}d-\p 2 n ,m) (30) 

n=— oo ' 

b(T,m)= 1 -^T £ / ^Hp 2 n )[rn + b(T, m )f (p 2 n )]d-\p 2 n , m ), (31) 

n=— oo ' 

where d{jp 2 w m) = p 1 A 2 {p 2 n ^ m) + u^C 2 ^^, m) + B 2 {j? w m). These equations are numerically 
solved. Note that the confining separable model quite facilitates summation over frequencies 
as well as the three-dimensional integration owing to the sufficient ultraviolet suppression 
guaranteed by the Gaussian functions fo{Pn) and fi(p 2 )- As a result, no renormalization is 
needed, which is the same as the zero-temperature case. 
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We now turn to the dressed scalar vertex. The finite temperature version of the inhomo- 
geneous ladder BSE for the dressed scalar vertex (Eq. (16)) with vanishing total momentum 
reads 

4 ±^ r (f 

r(p fc , 0; m) = 1 - -T J2 I J^9 2 D e JJ{p k - q n )^G{q n , m)T(q n , 0; m)G(q n , m) lv . (32) 

n=— oo ^ 1 

The solution of Eq. (32) takes the general form 

T(p n , 0, m) = F(p 2 n , m) ■ 1 + if • pG t (pl, m) + i-f 4 uj n G t (p 2 n , m). (33) 

The structure F{p 2 l) m) makes the leading-order contribution [30|] . In order to be consistent 
with the zero-temperature case and for simplicity as well, we shall keep only this term and 
expect it would represent well the temperature-dependent behavior of the dressed scalar 
vertex. Putting this ansatz T(p n ,0;m) = F(p^,m), the model gluon propagator (Eq. (26)) 
and the quark propagator (Eq. (28)) into Eq. (32), we find 

r> n ,0;m) = F(p 2 n ;m) = 1 + a(T, m)f (p 2 n ), (34) 

where a(T,m) satisfies the following equation 

16A)^ ^ M 3 <7^2^n, ,m w, 2m„ q 2 A 2 (q 2 ,m)+u; 2 n C 2 (qlm)-B 2 (qlm) 



a(T,m) = ^T £ J ^f (q 2 n )x[l+a(T,m)f (p 2 n )}x- 



[q 2 A 2 (ql, m) + u) 2 n C 2 {q 2 n , m) + B 2 (q 2 n , m)} 2 

(35) 

Having obtained the amplitudes A(q 2 ,m), B(q 2 ,m) and C(q 2 ,m) from Eqs. (29-31)), we 
can numerically solve a(T, m). The Gaussian function fo{Pn) appearing in the numerator of 
the integrand justifies an effective cutoff for summation over Matsubara frequencies as well 
as 3-D integration. 

Now, with the dressed quark propagator and dressed scalar vertex at finite temperature 
solved, the finite-temperature chiral susceptibility Eq. (8) can be finally expressed as 

^ f d 3 p 2 ~p 2 A 2 (p 2 n ,m) + uj 2 n C 2 {p 2 n) m)- B 2 (p 2 n ,m) 

(36) 

The factor "1" (corresponding to the bare part of the dressed scalar vertex) in the square 
bracket of the integrand would lead to quadratic divergence, which will be removed by 
subtracting the zero-temperature counterpart in our renormalization treatment. In order for 
the finite-temperature chiral susceptibility to approach a correct zero-temperature limit and 
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then the temperature-independent quadratic divergence to be canceled by the counterpart in 
the zero-temperature chiral susceptibility, one should employ the same numerical cut-off as 
that used in calculating x( m ) when calculating x(T, m), i.e. choose (^ 2 +o;^) max = Af, = A 2 ,, 
where At and A denote the numerical cut-off in calculating x(T, m) (Eq.( 36)) and x( m ) 
(Eq. (25)), respectively. 

C. Numerical results for renormalized chiral susceptibility 

Having obtained both the zero-temperature and finite-temperature chiral susceptibility, 
we are now in a position to calculate the renormalized chiral susceptibility. The numerical 
results are shown in Fig. 1 and Fig. 2, where Xr{T) are both scaled by T 2 and hence 
dimensionless. 




160 180 

T(MeV) 

FIG. 1: T-dependence of the renormalized chiral susceptibility in the chiral limit. 

From Fig. 1, one sees that in the chiral limit, the renormalized chiral susceptibility 
exhibits a very narrow, pronounced and in fact, divergent peak at the chiral critical tem- 



perature T c = 150 MeV, which is a typica 



driven by chiral symmetry restoration 101 ]. 



characteristic of phase transition of second order 



Fig. 2 shows a quite different picture at finite current quark masses (m = 0.0053 GeV 
corresponds to the physical value of the degenerate u, d-quark current mass in our model 



281]). There is no singular behavior any more: the peak of the susceptibility becomes 



not so sharp and pronounced as in the chiral limit; its height is greatly suppressed and 
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FIG. 2: T-dependence of the renormalized chiral susceptibility at different finite current quark 
masses, m = 0.0053 GeV corresponds to the physical quark mass in our model. 



evidently finite; and there is no unique temperature, but instead a range of finite width 
where the transition phenomenon takes place. Also notable are the decrease in the height 
of the susceptibility peak and that the susceptibility peak has a minor shift from 150 MeV 
to 155 MeV as the current quark mass increases. All of these observations agree well with 



lattice QCD results 



12j and support an analytic crossover involving a rapid change, 



as opposed to a jump, around the pseudo-critical temperature. 



IV. CONCLUSION AND DISCUSSION 



In the present work, we present a continuum investigation of the chiral susceptibility and 
try to say something about the character of the two-flavor QCD thermal transition from the 
temperature-dependent behavior of the susceptibility. We first derive a model-independent 
integral formula, which expresses the chiral susceptibility in terms of basic QFT objects: 
dressed propagator and vertex. After appropriately regularizing the additive quadratic di- 
vergence, we perform a model calculation of the renormalized chiral susceptibility within the 
DSE-BSE framework with a confining separable model gluon propagator which facilitates 
summation over Matsubara frequencies. Our model study shows that in the chiral limit, 
the renormalized chiral susceptibility exhibits a narrow, divergent peak at chiral critical 
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temperature T c = 150 MeV, signaling a second-order phase transition driven by chiral sym- 
metry restoration. On the other hand, at physical current quark masses, the renormalized 
chiral susceptibility features a crossover, which involves a rapid, but by no means singular, 
change in a temperature range of non-vanishing width. So, encouragingly, the lattice QCD 
observations 



121 ] are reproduced. 



Compared to other continuum model studies, like DSE model study [17j and NJL-type 



model study 15], the present work calculate the chiral susceptibility from an alternative 
definition (the former two model studies defined the susceptibility as the derivative of the 
generated mass with respect to the current quark mass and hence the quadratic divergence 
was circumvented) and in particular from a model-independent expression for the suscep- 
tibility. Within the present framework for the chiral susceptibility, the role of the dressed 
scalar vertex is consistently included into account. However, the present rainbow-ladder 
DSE-BSE model study still belongs to the class of mean-field approximations, like the two 
model studies mentioned above. It is interesting to investigate chiral susceptibility beyond 
mean-field approximations and this problem deserves further study. 
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